Electronic Transport in Disordered Bilayer and Trilayer Graphene 
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We present a detailed numerical study of the electronic transport properties of bilayer and trilayer 
graphene within a framework of single-electron tight-binding model. Various types of disorder 
are considered, such as resonant (hydrogen) impurities, vacancies, short- or long-range Gaussian 
random potentials, and Gaussian random nearest neighbor hopping. The algorithms are based on 
the numerical solution of the time-dependent Schrodinger equation and applied to calculate the 
density of states and conductivities (via the Kubo formula) of large samples containing millions of 
atoms. In the cases under consideration, far enough from the neutrality point, depending on the 
strength of disorders and the stacking sequence, a linear or sub linear electron-density dependent 
conductivity is found. The minimum conductivity a m i n ~ 2e 2 /h (per layer) at the charge neutrality 
point is the same for bilayer and trilayer graphene, independent of the type of the impurities, but 
the plateau of minimum conductivity around the neutrality point is only observed in the presence 
of resonant impurities or vacancies, originating from the formation of the impurity band. 

PACS numbers: 72.80.Vp, 73.22.Pr, 72.10.Fk 



I. INTRODUCTION 

Graphene is a subject of numerous investigations mo- 
tivated by its unique electronic and lattice properties, 
interesting both conceptually and for applications (for 
reviews, see Refs. llMlfil. Single layer graphene (SLG) is 
the two-dimensional crystalline form of carbon with a lin- 
ear electronic spectrum and chiral (A-B sublattice) sym- 
metry, whose extraordinary electron mobility and other 
unique features hold great promise for nanoscale electron- 
ics and photonics. Bilayer and trilayer graphenes, which 
are made out of two and three graphene planes, have also 
been produced by the mechanical friction and motivated 
a lot of researches on their transport properties 1 - 1 - - — . The 
charge-carrying quasiparticles in bilayer graphene (BLG) 
obey parabolic dispersion with non-zero mass, but retain 
a chiral nature similar to that in SLG (with the Berry 
phase 2-7T instead of 7T) 11 ' 12 . Furthermore, an electronic 
bandgap can be introduced in a dual gate BL G 15 ' 22 ' 41 " —, 
and it makes BLG very appealing from the point of view 
of applications. The trilayer graphene (TLG) is shown 
to have different electronic properties which is strongly 
dependent on the interlayer stacking sequenc e - 45 ' 46 . Nev- 
ertheless, graphene layers in real experiments always have 
different kinds of disorder, such as ripples, adatoms, ad- 
molecules, etc. One of the most important problems in 
the potential applications of graphene in electronics, is 
understanding the effect of these imperfections on the 
electronic structure and transport properties. 

The scattering theory for Dirac electrons in SLG is 
discussed in Refs. l47H5ll . Long-range scattering cen- 
ters are of special importance for transport properties 
of SLG, such as charge impurities 6 -! 5 ^— , ripples cre- 
ated long-range elastic deformations 7,56 , and resonant 
scattering centers 4 ^^^^^. Recently, the impact 
of charged impurity scattering on electronic transport 



in BLG have been investigated theoreticall y 17 ' 36 ' 37 and 
experimentally—. The linear density-dependent con- 
ductivity at high density and the minimum conduc- 
tivity behavior around the charge neutrality point are 
expecte d 17 ' 36 ' 37 and confirmed^ 8 ., but the experimental 
results also suggest that charged impurity scattering 
alone is not sufficient to explain the observed transport 
properties of pristine BLG on Si02 before potassium 
doping 38 . One possible explanation of the experimen- 
tal results might be the opening of a gap at the Dirac 
point in biased BLG 38 . On the other hand, some recent 
experimental 5 ^ and theoretical 6 ^— evidences appeared 
that the resonant scattering due to carbon-carbon bonds 
between organic admolecules and graphene (or by hydro- 
gen impurities which are almost equivalent to C-C bonds 
in a sense of electron scattering 6 *) is the main restricting 
factor for electron mobility in SLG on a substrate. These 
results suggest that the resonant impurity could also be 
the dominant factor of the transport properties of BLG 
and TLG. 

In the present paper, we study the effect of differ- 
ent types of impurities on the transport properties of 
graphene layers by direct numerical simulations in a 
framework of the single-electron tight-binding model. We 
consider four different types of defects: resonant ("hydro- 
gen") impurities, vacancies, Gaussian on-site potentials 
and Gaussian nearest carbon-carbon hoppings. The res- 
onant impurities/ vacancies and the centers of the Gaus- 
sian potentials/couplings are randomly introduced in the 
graphene layers. Our numerical calculations are based on 
the time-evolution metho d 65 ' 70 ' 71 , i.e., the time-evolution 
of the wave functions according to the Schrodinger equa- 
tion with additional averaging over a random superposi- 
tion of basis states. The main idea is that by performing 
Fourier transform of various correlation functions, such 
as the wave function-wave function and current-current 
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FIG. 1: Atomic structure of bilayer, ABA- and ABC-stacked trilayer graphene. 



correlation functions (Kubo formula), one can calculate 
the electronic structure and transport properties such as 
the density of states (DOS), quasieigenstates, ac (opti- 
cal) and dc conductivities. The details of the numerical 
method are presented in Ref. [65|. The advantages of the 
time-evolution method is that it allows us to carry out 
calculations for very large systems, up to hundreds of mil- 
lions of sites, with a computational effort that increases 
only linearly with the system size. 

The paper is organized as follows. Section II gives a de- 
scription of the tight-binding Hamiltonian of multilayer 
graphene. In section III, IV, V and VI, we focus on four 
different types of disorders respectively: resonant impuri- 
ties, vacancies, potential impurities, and nearest carbon- 
carbon hopping impurities. Finally a brief discussion is 
given in section VII. 



II. TIGHT-BINDING MODEL 

In general, the tight-binding Hamiltonian of multilayer 
graphene is given by 



N,, 



H = E E H 'u 
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where Hi is the Hamiltonian of SLG for Z'th layer and Hi 
describes the hopping between layers I and Z + 1. 
The single-layer Hamiltonian Hi is given by 



Hi = H + H v + H imp , 



(1) 



where Hq derives from the nearest neighbor hopping be- 
tween the carbon atoms: 



H = - 2 t 



(2) 



<%,j> 



H v denotes the on-site potential of the carbon atoms: 

H v = }^Vjcfci, (3) 



and Hi mp describes the resonant impurities (adatoms or 
admolecules) : 



H imp =£d^2 d+d. t + V ^2 ( d t c i + H - c - 



(4) 



The interlayer Hamiltonian H[ of bilayer graphene 
with AB Bernal stacking is given by 

H[ = -7i [ a t,j b l+hj + H - c - -73 [ b tj ai + 1 ^' + H - c - 

3 

(5) 

where a m i (b m .j) annihilates an electron on sublattice A 
(B), in plane m = I, I + 1, at site R (see the atomic struc- 
ture in Fig. [TJ. Thus, the second layer in BLG is rotated 
with respect to the first one by +120°. For the third 
layer there are two options: either the third carbon layer 
will be rotated with respect to the second layer by —120° 
(than it will be exactly under the first layer) or by +120°. 
In the first case we have ABA-stacked trilayer graphene, 
and in the second we have ABC-stacked (rhombohedral) 
graphene. The atomic structures of the ABA- and ABC- 
stacked trilayer graphene are shown in Fig. [TJ These 
stacked sequences can be extended to multilayers, i.e., 
the direct extension of ABA- and ABC-stacked sequences 
from trilayer to quartic-layer are ABAB and ABCD. The 
spin degree of freedom contributes only through a degen- 
eracy factor and, for simplicity, is omitted in Eq. ([TJ. 

The density of states is obtained by Fourier transfor- 
mation of the wave function at time zero and time t: 



1 f c 



(<p(0)\<p(t))dt, 



(G) 



where |y(0)) is an initial random superposition state 
of all the basis states and \<p(t)) = e~ lHt \(p (0)) is 
calculated numerically according to the time-dependent 
Schrodinger equation (we use units with h = 1). A de- 
tailed description of this method can be found in Refs. 
I65ll70l . The charge density is obtained by the intergral of 
the density of states, i.e., n e (E) = J Q E p (e) de. 

The static (dc) conductivity is calculated by using the 
Kubo formula 



" = - 1 A Tr {wSI' dt \ lJJ{t) + J{t)J] } 
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where J is the current operator and A is the sample area. 
The main idea of the calculation is to perform the time 
evolution of |<^(Q)). Then, we can extract not only the 
DOS but also the quasieigenstates \I> (e))^ 5 ., which are 
superpositions of degenerated energy-eigenstates. The 
conductivity at zero temperature can be represented as 



V 



dtRe [e- iet (ip (0)| Je lHt J \e)] , (8) 



where |e) is defined as 

\e) = 



1 



IM0)|*(e))l 



1*00) 



(9) 



The accuracy of the numerical results is mainly deter- 
mined by three factors: the time interval of the prop- 
agation, the total number of time steps, and the size 
of the sample. In the numerical calculations, the inte- 
grals in Eq. (O and (JSJ) are calculated using the Fast 
Fourier Transform (FFT) . According to the Nyquist sam- 
pling theorem, employing a sampling interval At = 
7r/maxi \Ei\, where Ei are the eigenenergies, is sufficient 
to cover the full range of eigenvalues. In practice, we do 
not know max^ | Ei | exactly but it is easy to compute an 
upperbound (for instance the 1-norm of H) such that At 
can be considered as fixed. 

In the present paper, the time evolution is calculated 
by the Chebyshev polynomial method, which has the 
same accuracy as the machine's precision independent of 
the value of time interval At. Alternatively, one could use 
Suzuki's product formula decomposition of the exponen- 
tial operators for the tight-binding Hamiltonian-iS, intro- 
ducing another time step that has to be (much) smaller 
than At to obtain accurate results^ .. In both cases, the 
accuracy of the energy eigenvalues is determined by the 
total number of the propagation time steps (N t ) that is 
the number of the data items used in the FFT. Eigenval- 
ues that differ less than AE = ir/N t At cannot be identi- 
fied properly. However, since AE is proportional to Nf 
we only have to extend the length of the calculation by a 
factor of two to increase the accuracy by the same factor. 

The third factor which determines the accuracy of our 
numerical results is the size of the sample. A sample 
with more sites in the real space will have more random 
coefficients in the initial state \ip(0)), providing a bet- 
ter statistical representation of the superposition of all 
energy eigenstates. This, however, is not a real issue in 
practice as it has be shown that the statistical fluctu- 
ations vanish with the inverse of the dimension of the 
Hilbcrt spaced, which for our problem, is proportional 
to the number of sites in the sample. A comparison of the 
DOS calculated from different samples size was shown in 
Ref. [65|, which clearly shows that larger sample size leads 
to better accuracy, and the result calculated from a SLG 
with 4096 x 4096 lattice sites matches very well with the 
analytical expression^ 5 -. More details on the numerical 
method itself can be found in Ref. [6^. The values of 
conductivity presented in this paper are normalized per 
layer and are expressed in units e 2 /h. 



Obviously, computer memory and CPU time evidently 
limit the size of the graphene system that can be simu- 
lated. The required CPU time is mainly determined by 
the number of operations to be performed on the state 
of the system, but this imposes no hard limit. However, 
the memory of the computer does. In the tight-binding 
approximation, a state \(p) of a sample consisting by N c 
atoms is represented by a complex- valued vector of length 
D = N c . For numerical accuracy (and in view of the 
large number of arithmetic operators performed), it is 
advisable to use 13—15 digit floating-point arithmetic 
(corresponding to 8 bytes per real number). Thus, to 
represent the state \ip) we need at least N c x 2 4 bytes. 
For example, for N c — 4096 x 4096 ~ 1.6 x 10 7 we need 
256 MB of memory to store a single arbitrary state \ip). 
This amount of memory is not a problem for the calcu- 
lation of DOS on a modest desktop PC or notebook, but 
it limits the calculation of the dc conductivity on such 
machines. To calculate one value of a (e) one needs stor- 
age of the corresponding quasieigenstate \e), and with 
typically 64 of such quasieigenstates in our simulations, 
a sample of N c = 4096 x 4096 sites requires at least 16 
GB memory for the storage, which is still reasonable for 
present-day computer equipment. 



III. RESONANT IMPURITIES 

Resonant impurities are introduced in reality by the 
formation of a chemical band between a carbon atom 
from graphene sheet and a carbon atom from an ad- 
sorbed organic molecule (CH 3 , C2H5, CH 2 OH), as well 
as H atoms 6 ^; vacancies are another option but in natu- 
ral graphene their concentration seems to be small. The 
adsorbates are described by the Hamiltonian Hi mp in 
Eq. ^j. From ab initio density functional theory (DFT) 
calculations^, it follows that the band parameters for 
various organic groups (and for hydrogen atoms) are al- 
most the same: V « 2t and w — i/16. The hybridiza- 
tion strength V being a factor 2 larger than t is in ac- 
cordance with the hybridization for hydrogen adatoms 
from Ref. 63, but the on-site energies are signifi- 
cantly smaller than the value Cd = 1.7eV used for H in 
Ref. |63j which makes our results for the transport prop- 
erties in SLG qualitatively different 64,65 . The adoption 
of these band parameters successfully explained the reso- 
nant scattering in SLG— and we continue to use them 
in the modeling of BLG and TLG. 

In Refs. I64l65l we used the algorithm presented 
in the previous section to calculate the dc conductiv- 
ity of SLG with resonant impurities or vacancies. We 
found that there is plateau of the order of the minimum 
conductivity 7 ^ 4e 2 /wh in the vicinity of the neutrality 
point, in agreement with theoretical expectations 7 ^. Be- 
yond the plateau, the conductivity is inversely propor- 
tional to the concentration of the impurities, and approxi- 
mately proportional to the carrier concentration n e . This 
is consistent with the approach based on the Boltzmann 
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FIG. 2: (Colour online) Top panel: DOS of bilayer graphene 
(71 = 73 = O.lt) with different concentrations of resonant 
impurities (ed = — 1/16, V = 2t) added on both layers. Mid- 
dle panel: Comparison of the conductivity of the BLG (line) 
and SLG (square) with the same concentration of resonant 
impurities. Bottom panel: Comparison of the conductivity 
of bilayer graphene (71 = 73 = O.lt) with the same amount 
of resonant impurities (ed = — 1/16, V = 2t) added on both 
layers (line, m) or only one layer (triangle, nu). Each layer in 
BLG contains 4096 x 4096 carbon atoms, and SLG contains 
6400 x 6400 carbon atoms. 



equation, which in the limit of resonant impurities with 
V -> 00, yields for the conductivityi&££i£&24 



(2e*/h) 



7r m 



e, 



D 



(10) 



carbon atom, and D is of order of the bandwidth. Note 
that for the case of the resonance shifted with respect to 
the neutrality point the consideration of Ref. leads to 
the dependence 



a oc (qo ± Uf hi fc^i?) 



(11) 



where ± corresponds to electron and hole doping, respec- 
tively, and R is the effective impurity radius. The Boltz- 
mann approach does not work near the neutrality point 
where quantum corrections are dominan t 57 > 72 i 73 . In the 
range of concentrations, where the Boltzmann approach 
is applicable, our numerical results of the conductivity of 
SLG as a function of energy fits very well to the depen- 
dence given by Eq. (fTTj^ ^. 

Electron scattering in BLG has been proven to dif- 
fer essentially from SLG in Ref. [T3 For a scattering 
potential with radius much smaller than the de Broglie 
wavelength of electrons, the phase shift of s-wave scat- 
tering Sq tends to a constant as k — >■ 0. Therefore, within 
the limit of applicability of the Boltzmann equation, the 
conductivity of a bilayer should be just linear in n e , in- 
stead of sublinear dependence (fTTjl for SLG. The differ- 
ence is that in SLG, due to vanishing DOS at the Dirac 
point, the scattering disappears at small wave vectors as 
So(k) oc 1/lnfci? (with In 2 kR on the order of 10 for typ- 
ical amounts of doping) for resonant and as So(k) oc kR 
for the nonresonant impurities. In contrast, in BLG there 
are no restrictions on the strength of the scattering and 
even the unitary limit (So = ir/2) can be reached at k = 0. 

However, these conclusions are based on the use of an 
approximate parabolic spectrum for the bilayer which is 
valid for the energy interval 



In the opposite case 



|£|«| 7 i|. 



|£|»|7i| 



(12) 



(13) 



where n e = E F /D 2 is the number of charge carriers per 



the effects of the interlayer hopping are negligible and 
one should expect a behavior of the conductivity similar 
to that of SLG. 

Our first set of numerical calculations of BLG are 
performed for similar concentrations of resonant impu- 
rities (m € [0.1%, 2%]) as those used for SLG in Refs. 
I64ll65l . The interlayer hopping parameters are taken as^ 
71 = 73 = O.K. As shown in Fig. [21 finite concentra- 
tions of the resonant impurities lead to the formation of 
a low energy impurity band (see increased DOS at low 
energies in Fig. [2]). The impurity band can host two 
electrons per impurity, and for impurity concentrations 
in the range of [0.1%, 2%], this leads to a plateau-shaped 
minimum of width 2rii in the conductivity vs. n e curves 
around the neutrality point. As one can see from the 
DOS in Fig. [5J even for n, = 0.1%, the width of the 
impurity band around the neutrality point is comparable 
to the limits of applicability of the parabolic approxima- 
tion for the spectrum (|12l) . therefore for the concentra- 
tions of the impurities presented in Fig. [5] one cannot 
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FIG. 3: (Colour online) Comparison of the DOS and conductivity of the SLG, BLG, TLG and QLG with the same concentration 
of resonant impurities (m = 0.5%, Sd = — i/16, V = 2t). 71 = 73 = O.lt in top panels and 71 = 0.5t, 73 = O.lt in bottom 
panels. SLG contains 6400 x 6400 carbon atoms, each layer in BLG, TLG and QLG contains 4096 x 4096, 3200 x 3200 and 
2400 x 2400 carbon atoms, respectively. 



use the theory^. For small electron concentrations we 
are beyond the limit of the Boltzmann theory at all, and 
for the larger electron concentration we are, rather, in 
the regime (fL3| so one could expect a sublinear behav- 
ior similar to that in SLG. Indeed, the conductivity of 
BLG as a function of charge density n e follows almost 
exactly the same dependence as for the SLG (see the di- 
rect comparisons of conductivities in Fig. [2|. That is, the 
density-dependence of conductivity in BLG is not linear 
but sublinear (ITT]) as in SLG. Actually, as shown in Fig. 
[21 the sublinear dependence is quite general for multilayer 
graphene, i.e., it is also true for trilayer and quartic-layer 
graphene with the same concentration of resonant impu- 
rities, independent on the stacking sequence, which is, of 
course, not surprising assuming that the condition (|T5|) 
holds. Here for the trilayer (quartic-layer) we consider 
two types of stacking sequence: ABA (ABAB) and ABC 
(ABCD). This general property of the conductivity can 
be easily understood by comparison of their DOS in Fig. 
[3l The DOS of single-layer, bilayer, trilayer and quartic- 
layer graphene are exactly the same except near the edge 
of the spectrum, indicating the similar band structure, 
independent on the number of layers and stacking se- 



quence. In fact, since the couplings between the carbon 
atoms and organic admolecules are twenty times larger 
than the interlayer coupling (V — 2O71) in our model, 
the unique bonds generated by the relevant weaker inter- 
layer interactions are more easily to be destroyed by the 
impurity bonds generated by the much stronger adsorbed 
resonant impurities. 

In order to check the symmetry of the presence of 
the impurities, we limit the adsorption of organic ad- 
molecules to one layer of BLG (nu, case II). To com- 
pare the results of the adsorption on both sides (n.j , case 
I), we fix the total number of resonant impurities and 
therefore the concentration on one layer (case II) is dou- 
bled (nij = N imp /N carbon _ m _ one j ayer = 2ni). We find 
(see last panel in Fig. [2]) that for the low concentration 
(rij < 0.5%), the electron-density dependence of the con- 
ductivity in BLG follows the same law in both cases; at 
high concentration (rii > 1%), the conductivity in case 
II is larger than in case I. This is because in case II the 
difference of mobility of electron in the two layers, with 
or without impurities, is larger than in the case of large 
concentrations of adsorbed admolecules. 

Next we consider the region of parameters which can 
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FIG. 4: (Colour online) Comparison of the DOS of bilayer 
graphene with different interlayer interactions: 71 = O.lt, 
0.2i, and 0.5t (73 is fixed as O.lt). Inner panel: normalized 
DOS and energy in units of I/71 and 71. Each layer in BLG 
contains 4096 x 4096 carbon atoms. 



be described by the Boltzmann equation plus parabolic 
spectrum^. In BLG, the approximations of massive 
valence and conduction bands with zero gap: E {k) = 
±H 2 k 2 /2m*, where the effective mass is given by m* = 
ji/2Vp, are only true in the low-energy dispersion close 
to the neutrality point. There are two ways to place the 
impurity bands within the region of low-energy dispersion 
(|12p : decreasing the concentration m of impurities or ex- 
panding the quadratic band by increasing 71. Smaller 
concentration of impurities leads to less random states 
for the averaging in the Kubo formula of Eq. ([5]) , which 
means that it is numerically more expensive because we 
need to extend the sample size to keep the same accuracy. 
Therefore increasing 71 is computationally more conve- 
nient from the point of view of CPU time and physical 
memory; one can assume that physical results should be 
the same: it is only the ratio rii/71 which is important. 

In Fig. |4j we compare numerical results of DOS of BLG 
with different band parameter 71: O.lt, 0.2t and 0.5t (73 is 
fixed as O.lt). One can see that the width of the parabolic 
band with the energy-independent local density of states 
proportional to 71, and the normalized energy (in units 
of 71) dependencies of DOS (in the units of I/71) within 
the parabolic band are consistent for different 71 (sec the 
inner panel of Fig. 0|) . Therefore we can simply use 71 
with the value of 0.5t instead of O.lt to extend the width 
of the parabolic band approximation without changing 
the structure of the spectrum. The numerical results 
form a system with m times larger of 71 , are qualitatively 
comparable to those for a system of 1/m times smaller 
concentration rii of impurities. 

The numerical results for the DOS and conductivities 
of BLG and TLG in the presence of resonant impurity 
with larger interlayer interactions (71 = 0.5i,) are shown 
in Fig. [3l We see that for an impurity concentration 
of rii = 0.5%, the impurity band is located around the 



neutral point and far from the edge of the quadratic 
band {\E\ < 0.5t). In the region of the impurity band 
(\n e \ < rii — 0.5%), there is a plateau in the order of 
2e 2 /h (per layer) in BLG, as well as in TLG. This values 
is slightly larger than the minimum conductivity 4e 2 /nh 
of SLG. It is worthwhile to note that an explanation of 
the origin of plateau around the neutrality point is be- 
yond the applicability of Boltzmann equation, just as in 
the case of SL G 64 i 65 . Analyzing experimental data of the 
plateau width (similar to the analysis for N2O4 acceptor 
states in Ref. \7a ) can therefore yield an independent es- 
timate of the impurity concentration, both in single-layer 
and multilayer graphene. Within the parabolic band but 
beyond the impurity band, the conductivities in BLG 
and ABA-stacked TLG exhibit very well the linear de- 
pendence on the charge density n e . The ABC-stacked 
TLG is different from the others because of its unique 
band structures with a cubic touching of the bands^ (see 
the difference of DOS in Fig. [3]). 

Finally we check the role of 73 in the conductivity 
of BLG. Theoretically, the influence of 73 to the band 
structure is negligible, and so it is for the conductivity. 
This is confirmed by our numerical results in Fig. [5] 
For the fixed concentration of impurities rii — 0.5% with 
71 = 0M, the values of the conductivity corresponding 
to the same electron concentration n e are quite close for 
73 = 0, O.li, and 0M. 



IV. VACANCIES 

A vacancy in a graphene sheet can be regarded as an 
atom (lattice point) with an on-site energy v — > 00 or 
with its hopping parameters to other sites being zero. 
In the numerical simulation, the simplest way to imple- 
ment a vacancy is to remove the atom at the vacancy 
site. Introducing vacancies in SLG will create a zero en- 
ergy modes (midgap state) 48 ' 61 i 62 i 64 ' 65 . The exact ana- 
lytical wave function associated with the zero mode in- 
duced by a single vacancy in SLG was obtained in Ref. 
l6ll showing a quasilocalized character with the ampli- 
tude of the wave function decaying as inverse distance to 
the vacancy. SLG with a finite concentration of vacan- 
cies was studied numerically in Refs. H62.64.65,63-[6l 
It was shown that the number of the midgap states in- 
creases with the concentration of the vacancie o 48 ' 62 i 64 ' 65 , 
and quasieigenstates are also quasilocalized around the 
vacancies^. The inclusion of vacancies brings an increase 
of spectral weight to the surrounding of the Dirac point 
(E = 0) and smears the van Hove singularitie a 48 ' 62 i 64 ' 65 . 
The effect of the vacancies on the transport properties 
of SLG is quite similar to that of the adsorbed organic 
molecules. The main difference is the position of the im- 
purity band in the spectrum: its center is located at the 
neutrality point in the presence of vacancies, whereas it is 
biased in the presence of realistic resonant impurities be- 
cause of the nonzero on-site potential on the organic car- 
bon (or hydrogen) atom. The vacancy band contributes 
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FIG. 5: (Colour online) DOS and conductivity of bilayer graphene with resonant impurities (Ed ~ — 1/16, V — 2t) added on 
both layers. The interlayer parameter 71 is fixed as 0.5t, and 73 = (black line), O.li (red dashed line) and 0.5t (green dot 
line). Each layer contains 4096 x 4096 carbon atoms. 




FIG. 6: (Colour online) Comparison of the DOS and conductivity of the SLG, BLG and TLG with the same concentration of 
vacancies (n x = 0.5%). The parameters of the interlayer coupling are 71 = 0.5t and 73 = O.lt. SLG contains 6400 x 6400 carbon 
atoms, each layer in BLG and TLG contains 4096 x 4096 and 3200 x 3200 carbon atoms (sites), respectively. 



to the conductivity and leads to a plateau of minimum 
conductivity in the midgap region. The width of the 
plateau is 2n x (n x is the concentration of the vacancies) 
in the conductivity vs. n e curves around the neutrality 
point, showing the same dependence (2m) as the case 
of resonant impuritie o 64 i 65 . For the range of concentra- 
tions where the Boltzmann approach is applicable, the 
conductivity of SLG as a function of energy fits very well 
to the dependence given by Eq. (ITU , with q = for the 
vacancies and qo 7^ for the resonant impuritie c 64 ' 65 . 

Previous studies of the vacancies in BLG focused 
mainly on the properties of the local density of states 
(LDOS) around a single or a pair of vacancies, and it 
was shown that the LDOS in the neighboring lattice sites 
of the impurity site is normally enhanced, depending on 
the lattice site (A or B sublattices) of the vacanc y 76 ' 77 . 
Recently, a new type of zeromodel state in BLG is found 
in Ref. [78|, in the absence of a gap it is quasilocalized 



in one of the layers and delocalized in the other, and in 
the presence of a gap it becomes fully localized inside the 
gap. These observations are different from SLG, where 
the impurity state is insensitive to the position of vacan- 
cies. The differences in the spectrum of LDOS around 
the vacancies in SLG and BLG lead to different electron- 
density (Fermi energy) dependence of the conductivity. 
As the vacancies and resonant impurities have similar ef- 
fects on the electronic structure and transport properties 
in SL G 64 ' 65 , it suggests that their contributions to the 
bilayer and trilayer graphene should also be comparable. 
We consider here the results for the vacancies in the range 
that the Boltzmann approach is applicable. In Fig. [5J we 
show the numerical results of the DOS and conductivi- 
ties of SLG, BLG and TLG with fixed concentration of 
vacancies (n x = 0.5%). The parameters of the interlayer 
coupling are 71 = 0.5i and 73 = O.lt. These results are 
directly comparable with the results of the same concen- 
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FIG. 7: (Colour online) Contour plot of the on-site potentials in the central part of a graphene layer (4096 x 4096) with 
short-range (A = 3t, d = 0.65a, P v — 0.5%) or long-range (A = It, d = 5a, P v — 0.1%) Gaussian potential. 




FIG. 8: (Colour online) DOS and conductivity of bilayer graphene (71 = 73 = O.lt) with short-range (A = 3t, d = 0.65a) or 
long-range (A = It, d = 5a) Gaussian potential. Each layer contains 4096 x 4096 carbon atoms. 



tration of resonant impurities represented in Fig. |3l and resonant impurities, 
demonstrate similar density-dependence of the conduc- 
tivities, just as we expected. For conciseness, we do not 
discuss these vacancies as their effect on the transport 
properties of graphene are quite similar to those of the 
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FIG. 9: (Colour online) DOS and conductivity of bilayer graphene (71 = 0.5t, 73 = O.li) with long-range (A = It, d = 5a) or 
short-range (A = St, d = 0.65a) Gaussian potential. Each layer contains 4096 x 4096 carbon atoms. 



V. GAUSSIAN POTENTIAL 



The impurities in the Hamiltonian of Eq. (JTJ are rep- 
resented by random on-site potentials. Short-rang and 
long-range Gaussian potentials are given by 



NT, 



= E u * 



exp 



k=l 



|£i ~ £fc| 

2d 2 



(14) 



where Nf mp is the number of the Gaussian centers, which 
are chosen randomly distributed on the carbon atoms, 
Uk is uniformly random in the range [—A, A] and d is 
interpreted as the effective potential radius. The typical 
values of d used in our model are d = 0.65a and 5a for 
short- and long-range Gaussian potential, respectively. 
Here a is the carbon-carbon distance in the monolayer 
graphene. The value of Nf mp is characterized by the value 
P v — Nj mp /N, where N is the total number of carbon 
atoms of the sample. A typic contour plot of the on-site 
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FIG. 10: (Colour online) DOS and conductivity of bilayer graphene (71 = 73 = O.lt) with short-range (At = t, dt = 0.65a) or 
long-range (At = O.bt, d t = 5a) Gaussian hopping. Each layer contains 4096 x 4096 carbon atoms. 



potentials in the central part of a graphene layer with 
short- or long-range Gaussian potential is shown in Fig. 
UJ The sum in Eq. (IT4l is limited to the sites in the 
same layer, i.e., we do not consider the overlapping of 
the Gaussian distribution in different layers. 

Numerical results of the density of states and dc con- 
ductivities of BLG (71 = 73 = O.lt) with short- (A = M, 
d = 0.65a) and long-range (A = It, d = 5a) Gaussian 
potentials are shown in Fig. [51 Similar to the case of 
resonant impurities, the singularities in the spectrum are 
also suppressed in the presence of random potentials, and 
the conductivity as a function of charge density follows 
a sublinear dependence. The difference is that there is 
no impurity band around the neutrality point (see the 
DOS in Fig. [5J. This leads to totally different trans- 
port properties: no plateau around the Dirac point in 
the conductivity vs. n e curves. 

Similar to the case of resonant impurities, the regime 
of parabolic band in BLG expands by increasing 71 from 
O.lt to 0.5t, the results being shown in Fig. [HI Now 
the difference of transport properties in BLG with short- 
and long-range Gaussian potentials are more significant 
within the parabolic band: the density-dependence of 
conductivity is sublinear in the case of short-range, but 



linear in the case of long-range potentials. Actually, these 
sublinear and linear dependencies are also observed in 
TLG, independent on the stacking sequence (see Fig. [§]). 

The same value of the minimum conductivity (cr m i n « 
2e 2 /h) at the charge neutrality point is observed for both 
BLG and TLG with 71 = 0.5i. As we discussed in the 
case of resonant impurities, the adoption of larger 71 
is equivalent to the use of smaller disorder, and there- 
fore our results indicate that the minimum conductivity 
in order of a m i n ~ 2e 2 /h is common in BLG and TLG 
with small concentration of random Gaussian potentials. 
These numerical results are consistent with the analytical 
result for BLG in Ref. E3- 



VI. GAUSSIAN HOPPING 

The origin of disorder in the nearest neighbor coupling 
could be substitutional impurities like N or B instead of 
C, or distortions of graphene sheet. To be specific, we 
introduce the disorder in the hopping by a Gaussian dis- 
tribution in a similar way as random Gaussian potential, 
namely, the distribution of the nearest neighbor hopping 
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FIG. 11: (Colour online) DOS and conductivity of bilayer graphene (71 = 73 = O.lt) with short-range (A t = t, d t — 0.65a) or 
long-range (At = O.bt, d t = 5a) Gaussian hopping. Each layer contains 4096 x 4096 carbon atoms. 



parameters reads 

Ui =t+i: T k exp (^- Ir ' + r ^ 2rfc|2 ) , (15) 

where N* is the number of the Gaussian centers, Tk 
is uniformly random in the range [—At, A t ] and d t is in- 
terpreted as the effective screening length. Similarly, the 
typical values of dt are the same as for the Gaussian po- 
tential, i.e., dt = 0.65a and 5a for short- and long-range 



Gaussian random hopping, respectively, and the values of 
N- mp are characterized by the value P t = N- mp /N. Sim- 
ilar as in Eq. (fH|) . the sum in Eq. (fT5")) does not include 
the overlapping of the Gaussian distribution in different 
layers. 

Like in the case of Gaussian potentials, the presence 
of random Gaussian hopping in BLG and TLG also sup- 
presses the Van Hove singularities in the spectrum, but 
does not introduce a new impurity band (midgap states) 
and there is also no plateau in the conductivity vs. elec- 
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tron density curves (see Fig. [TUland [TT|) . The unique fea- 
ture characteristic for the presence of random Gaussian 
hopping is that in the region near the neutrality point, 
the conductivity is always linearly dependent on the elec- 
tron density, with no influence from the concentration 
of Gaussian centers (different Pt in Fig. ITU1) . range of 
Gaussian coupling (dt = 0.65a or 5a), strength of the in- 
terlayer coupling (71 = O.lt in Fig. ITTJ1 and 0.5t in Fig. 
Ill[) . number of layers (bilayer or trilayer) and stacking se- 
quence (ABA or ABC in TLG). The differences of short- 
or long-range cases are only obvious in the energy re- 
gion far from the neutrality point (high concentration of 
charge density) : the increase of conductivity as a function 
of charge density is monotonic only for the long-range dis- 
order. Furthermore, like in the case of random Gaussian 
potential, a common minimum conductivity in the order 
of 2e 2 /h on charge neutrality point is also observed for 
both BLG and TLG. 



VII. DISCUSSION AND CONCLUSIONS 

We have presented a detailed numerical study of the 
electronic transport properties of bilayer and trilayer 
graphene within the framework of a noninteracting tight- 
binding model. Various realistic types of disorder are 
considered, such as resonant impurities, vacancies, ran- 
dom Gaussian on-site potentials, and random Gaussian 
hopping between nearest carbon atoms. Our results 
give a consistent picture of the electronic structure and 
transport properties of bilayer and trilayer graphene in 
a broad range of concentration of impurities or other 
sources of disorder. Linear or sublinear electron-density 
dependent conductivity at high enough density is ob- 
served, depending on the type and strength of the disor- 
der and the stacking sequence. The minimum conductiv- 
ity cr m i n ?» 2e 2 /h (per layer) on charge neutrality point is 
common for BLG and TLG, independent of the type of 
the impurities, but the plateau of minimum conductiv- 
ity around the neutrality point is unique when resonant 
impurities or vacancies are present. 

In the presence of resonant impurities or vacancies, the 
dependence of the conductivity as a function of electron 
density is affected by the relevant width of the impurity 
band and the band created by the interlayer hopping. 
Using BLG with vacancies as an example: introducing 



n p = n e (71) = J Q 71 p (e) de as the density of electrons on 
the boundary of the parabolic band, and considering the 
case that the concentration of vacancies (n x ) is smaller 
than rip, i.e., the impurity band is within the region of 
the parabolic conduction band, there are three regions of 
electron-density dependence of the conductivity: 

(i) |n e | < n x , a central minimum conductivity plateau 
(2e 2 /h per layer) with width equals to 2n x ; 

(ii) n x < \n e \ < n p , linear dependence, as predicted by 
the analytical treatment using the Boltzmann equation 
for parabolic spectrum^; 

(hi) n x > n pi sublinear dependence, as the effects of 
the interlayer hopping are negligible in this region and 
one should expect a behavior of the conductivity similar 
to that of SLG. 

On the opposite case n x > n p , region (ii) simply dis- 
appears and therefore we can only observe the minimum 
conductivity plateau and sublinear dependence on the 
high concentration of electron densities. Actually, the 
sublinear dependence beyond the parabolic band is a gen- 
eral property of SLG, BLG and MLG with large enough 
concentration of resonant impurities or vacancies, inde- 
pendent on the number of layers and the stacking se- 
quence. 

In the presence of random Gaussian on-site potentials, 
the electron-density dependences of conductivity of BLG 
or TLG arc sublinear and linear in the low concentration 
of charges, for short- and long-range disorders, respec- 
tively but are always sublinear in the high concentra- 
tion. On the other hand, in the case of random Gaus- 
sian carbon-carbon couplings, the density-dependence of 
conductivity in the region close to the neutrality point is 
more simple: there is only a linear dependence, with no 
effect of the strength and range of disorder, the number 
of layers and stacking sequence. 

Note added: After this paper was submitted, a paper 
which also discusses the effect of resonant scatterers on 
the dc conductivity of single-layer and bilayer graphene 
appeared^, with results that are consistent with ours. 
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